/*
 * level1.c
 *
 * Copyright (c) 2009, Eugen Stoian <stoian.e@gmail.com>
 * 
 * This library is free software; you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the Free
 * Software Foundation; either version 2.1 of the License, or (at your option)
 * any later version.
 * 
 * This library is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
 * details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
 * http://www.fsf.org/licensing/licenses/lgpl.txt
 */
#include "level1.h"
#include <math.h>
/*
 * SUM
 */
void ssum(int n, const float* x, const float* y, float* result){
	int i;
	for(i = 0; i < n; ++i)
		result[i] = x[i]+y[i];
}
void dsum(int n, const double* x, const double* y, double* result){
	int i;
	for(i = 0; i < n; ++i)
		result[i] = x[i]+y[i];

}/*
 * ASUM
 */
void sasum(int n, const float* x, const float* y, float* result){
	int i;
	for(i = 0; i < n; ++i)
		result[i] = fabs(x[i])+fabs(y[i]);

}
void dasum(int n, const double* x, const double* y, double* result){
	int i;
	for(i = 0; i < n; ++i)
		result[i] = abs(x[i])+abs(y[i]);
}
/*
 * DOT PRODUCT
 */
float sdot(int n, const float* x, const float* y){
	float result = 0.0f;
	int i;
	for(i = 0; i < n; ++i){
		result += x[i]*y[i];
	}
	return result;
}
double ddot(int n, const double* x, const double* y){
	double result = 0.0;
	int i;
	for(i = 0; i < n; ++i){
		result += x[i]*y[i];
	}
	return result;
}

/*
 * EUCLIDEAN NORM
 */
float snrm2(int n, const float* x){
	return sqrtf(sdot(n,x,x));
}
double dnrm2(int n, const double* x){
	return sqrt(ddot(n,x,x));
}

/*
 * y <- a*x+y
 */
void saxpy(int n, float a, const float* x, float* y){
	int i;
	for(i=0;i<n;++i){
		y[i]+=a*x[i];
	}
}
void daxpy(int n, float a, const double* x, double* y){
	int i;
	for(i=0;i<n;++i){
		y[i]+=a*x[i];
	}

}

/*
 * x <- a*x
 */
void sscale(int n, float a, float* x){
	int i;
	for(i=0;i<n;++i){
		x[i]*=a;
	}
}
void dscale(int n, double a, double* x){
	int i;
	for(i=0;i<n;++i){
		x[i]*=a;
	}
}

/*
 * IAMAX
 */
int samax(int n, const float* x){
	float amax = 0.0f;
	int imax = 0;
	int i;
	for(i = 0; i < n; ++i)
		if(amax<fabs(x[i])){
			amax = fabs(x[i]);
			imax = i;
		}
	return i;
}

int damax(int n, const double* x){

	double amax = 0.0;
	int imax = 0;
	int i;
	for(i = 0; i < n; ++i)
		if(amax<abs(x[i])){
			amax = abs(x[i]);
			imax = i;
		}
	return i;
}
